Introduction

library(tidyverse)
library(raster)
library(sf)
library(enerscape)

#' Function to plot the areas
#' 
#' @param x is a raster
#' @param poly is a shapefile polygon
plot_area <- function(x, poly, mask = NULL, col = NULL, void = FALSE) {
  if (is.null(col)) {
    plot(x, axes = !void, box = !void)
  } else {
    plot(x, col = col, axes = !void, box = !void)
  }
  if (!is.null(mask)) {
    plot(mask, add = T, legend = FALSE,
         col = adjustcolor("white", alpha.f = 0.5), 
         axes = !void, 
         box = !void)
  }
  plot(poly$geometry, add = TRUE)
}

#' ARC model predictions
#'
#' @param mass body mass of the animal (kg)
#' @param slope incline of the terrain (degree)
ARC <- function(mass, slope) {
  E_ar <- 8 * mass^(-0.34)
  E_mec <- 50 * (1 + sin((2 * slope - 74) / 180 * pi)) * mass^(-0.12)
  E <- E_ar + E_mec
  return(E)
}

pontzer %>% 
  as_tibble %>% 
  mutate(ARC = ARC(Mass, Incline)) %>% 
  ggplot() +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  geom_point(aes(ARC, Cost.of.Transport, col = Incline)) +
  scale_color_gradient(low = "steelblue", high = "tomato") +
  scale_x_log10(n.breaks = 11) +
  scale_y_log10(n.breaks = 11) +
  xlab("ARC predicitons") +
  ylab(expression("E"["COT"]*" (J kg"^"-1"*"m"^"-1"*")")) +
  theme_bw()

dem <- raster("../data/dem-vignette.tif")
res(dem)

## [1] 100 100

crs(dem)

## CRS arguments:
##  +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

en <- enerscape(dem, 50)

##  - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)

##   | Calculating slope

##   | Calculating work

##   | Calculating conductance (1 / work)

##  - Do not use slope with negative values for distance calculations

extr <- en_extrapolation(en, plot = FALSE)
extr

## $`Mass extrapolated`
## [1] FALSE
## 
## $`Slope extrapolated`
## [1] TRUE
## 
## $`Slope extrapolation`
## class      : RasterLayer 
## dimensions : 367, 509, 186803  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : Slope 
## values     : 1, 1  (min, max)

plot(dem)
plot(extr$`Slope extrapolation`, 
     add = TRUE, 
     col = adjustcolor("blue", alpha.f = 0.3),
     legend = FALSE)
arrows(888523, 4674275, 878791, 4671549)

lcp <- enerscape::en_lcp(en, simulate_random_points = TRUE)

## 1 of 10 random points ...

## 2 of 10 random points ...

## 3 of 10 random points ...

## 4 of 10 random points ...

## 5 of 10 random points ...

## 6 of 10 random points ...

## 7 of 10 random points ...

## 8 of 10 random points ...

## 9 of 10 random points ...

## 10 of 10 random points ...

walk <- enerscape::en_passage(en,
                              simulate_random_points = TRUE,
                              theta = 4)

## 1 of 10random points ...

## 2 of 10random points ...

## 3 of 10random points ...

## 4 of 10random points ...

## 5 of 10random points ...

## 6 of 10random points ...

## 7 of 10random points ...

## 8 of 10random points ...

## 9 of 10random points ...

## 10 of 10random points ...

walk

## $Origins
##            x       y
##  [1,] 881090 4694830
##  [2,] 883990 4676830
##  [3,] 870990 4679430
##  [4,] 857490 4690330
##  [5,] 897890 4688530
##  [6,] 876390 4692730
##  [7,] 880990 4660130
##  [8,] 882590 4692530
##  [9,] 886190 4671730
## [10,] 868990 4679330
## 
## $Destinations
##            x       y
##  [1,] 854590 4672030
##  [2,] 874290 4679530
##  [3,] 874990 4692030
##  [4,] 880490 4677330
##  [5,] 883590 4669330
##  [6,] 867490 4661530
##  [7,] 893390 4676730
##  [8,] 868090 4689830
##  [9,] 893390 4671230
## [10,] 869490 4662730
## 
## $Passages
## class      : RasterStack 
## dimensions : 367, 509, 186803, 10  (nrow, ncol, ncell, nlayers)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## names      : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10 
## min values :       0,       0,       0,       0,       0,       0,       0,       0,       0,        0 
## max values :       0,       0,       0,       0,       0,       0,       0,       0,       0,        0 
## 
## 
## $`Cumulative net passage`
## class      : RasterLayer 
## dimensions : 367, 509, 186803  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : NA, NA  (min, max)

Landscape connectivity for the Marsican bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park.

Here, I assess the overall landscape connectivity for the Mariscan brown bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park (SVRP) in central Italy parcosirentevelino.it. In particular, I derive the energy landscape for the area and use it as a “resistance matrix” for animal movement, i.e. assuming that the animal would move across low energy landscape trajectories in order to minimize energy costs.

The required input for enerscape using the ARC model are:

• a digital elevation model (DEM) of the area of interest (m)
• body mass of the animal (kg)

I load the DEM and the park polygon data:

dem <- raster("../data/dem-vignette.tif")
park <- st_read("../data/SVRP-polygon.shp")

## Reading layer `SVRP-polygon' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/SVRP-polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 30 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 857140.5 ymin: 4665004 xmax: 897975.7 ymax: 4691681
## CRS:            unknown

m <- mask(dem, park, inverse = TRUE)
plot_area(dem, park, m, void = TRUE)
title("Elevation (m)")
scalebar(10000, 
         type = "bar", 
         divs = 2, 
         label = as.character(c(0, 5, 10)),
         below = "kilometers")

I calculate the energy landscape for a bear of 140 kg:

en <- enerscape(dem, 140, "kcal")

##  - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)

##   | Calculating slope

##   | Calculating work

##   | Calculating conductance (1 / work)

##  - Do not use slope with negative values for distance calculations

plot(en$rasters, axes = FALSE, box = FALSE)

Or plotting only the cost of travel per cell:

plot_area(en$rasters$Work, park, m, col = topo.colors(100), void = TRUE)
title("Energy cost of travel (kcal)")

I choose two points for which to compute the least-cost path and the connectivity using Circuitscape:

p <- data.frame(x = c(877367, 882653),
                y = c(4674192, 4677413))
plot_area(dem, park, m, void = TRUE)
points(p, pch = 20)

lcp <- en_lcp(en, p[1, ], p[2, ])

circuitscape_skeleton(en,
                      file.path("/home",
                                "eb97ziwi",
                                "Proj",
                                "enerscape-paper",
                                "data"),
                      p)

And I run Circuitscape in Julia.

julia> using Omniscape
julia> run_omniscape("/home/eb97ziwi/Proj/enerscape-paper/data/circuitscape.ini")

## [ Info: 2020-12-19 08:31:24 : Precision used: Double
## [ Info: 2020-12-19 08:31:27 : Reading maps
## [ Info: 2020-12-19 08:31:28 : Resistance/Conductance map has 186803 nodes
## [ Info: 2020-12-19 08:31:36 : Solver used: AMG accelerated by CG
## [ Info: 2020-12-19 08:31:36 : Graph has 186803 nodes, 2 focal points and 1 connected components
## [ Info: 2020-12-19 08:31:36 : Total number of pair solves = 1
## [ Info: 2020-12-19 08:31:37 : Time taken to construct preconditioner = 1.683233535 seconds
## [ Info: 2020-12-19 08:31:37 : Time taken to construct local nodemap = 0.008903604 seconds
## [ Info: 2020-12-19 08:31:37 : Solving pair 1 of 1
## [ Info: 2020-12-19 08:31:38 : Time taken to solve linear system = 0.526888286 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to write current maps = 0.559248897 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to complete job = 15.111858446
## 3×3 Array{Float64,2}:
## 0.0   1.0      2.0
## 1.0   0.0     31.8625
## 2.0  31.8625   0.0

Sometimes, few cells with extremely high current/connectivity can make the output map hard to read, as most of the connectivity values are squeezed into a narrow range. To avoid this, the value of extremely high outliers can be forced to an upper limit. For example, I force outliers to the 0.999 quantile value of the connectivity and standardize the Circuitscape output between 0 and 1.

cs <- raster("../data/circuitscape_cum_curmap.asc")
cs <- log10(cs)
cs <- cs - min(values(cs), na.rm = TRUE)
cs[cs >= quantile(cs, 0.999)] <- quantile(cs, 0.999) #remove far outliers
cs <- cs / max(values(cs), na.rm = TRUE)
plot_area(cs, 
          park, 
          m, 
          col = viridis::inferno(100), 
          void = TRUE)
points(p, pch = 3)
lines(lcp$Paths)
title("Connectivity between locations")

To compute the overall landscape connectivity, I use the Omniscape algorithm, which runs Circuitscape iteratively on a moving window. I initialize the omniscape.ini file to be passed to Julia:

omniscape_skeleton(en,
                   file.path("/home",
                             "eb97ziwi",
                             "Proj",
                             "enerscape-paper",
                             "data"),
                   radius = 10)

And I run the Omniscape algorithm in a Julia terminal. Don’t forget to specify the number of parallel threads: https://docs.circuitscape.org/Omniscape.jl/stable/usage/#Parallel-Processing.

julia> using Omniscape
julia> run_omniscape("/home/eb97ziwi/Proj/enerscape-paper/data/omniscape.ini")

## Starting up Omniscape. Using 7 workers in parallel. Using double precision...
## Solving moving window targets...
## Progress: 100% | Time: 0:03:14
## Done!
## Time taken to complete job: 196.3408 seconds
## Outputs written to
## /home/eb97ziwi//home/eb97ziwi/Proj/enerscape-paper/data/omniscape_2
## (Union{Missing, Float64}[0.5152914765395912 0.5884433266779532 … 0.719651357069539 0.8600629588479131;
## 0.4834003129933446 0.6594116484239914 … 0.635221929581483 0.7110678186022065;
## … ;
## 0.2127925084620419 0.2067469779201847 … 0.39577453680538544 0.46681040575465904;
## 0.17382917047741175 0.2093726458198999 … 0.31763823508216504 0.3738227777731409]
## Union{Missing, Float64}[1.0885252650071 1.136985249629082 … 1.169180462521136 1.1104937306261533;
## 1.053217843455199 1.0090282247541849 … 0.9301240063409788 1.1859951520336032;
## … ;
## 1.1798282079383584 0.9871029439774996 … 0.9873244151514602 1.2432973368603717;
## 1.1120398150491893 1.168088759875787 … 1.0583987239714665 1.0877633126993709])

I remove outliers and standardize the Omniscape output between 0 and 1.

os <- raster("../data/omniscape_2/cum_currmap.tif")
os <- os - min(values(os), na.rm = TRUE)
os[os >= quantile(os, 0.999)] <- quantile(os, 0.999) #remove far outliers
os <- os / max(values(os), na.rm = TRUE)
plot_area(os, 
          park, 
          m, 
          col = viridis::inferno(100), 
          void = TRUE)
title("Landscape connectivity")

Mols Bjerge

• Tree cover density: Copernicus
• NDVI: Sentinel 2
• DEM: kortforsyningen

Energy costs of cycling L’Eroica

Here, I study how body weight (kg) and cycling speed (km/h) influence the energy required by participants to cycle L’Eroica event. In particular, I derive the energy landscape for the route and sum the energy costs of all cells.

I load the DEM and the route polygon data:

dem <- raster("../data/eroica/eroica.tif")
route <- st_read("../data/eroica/route.shp") %>% 
  st_transform(crs(dem))

## Reading layer `route' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/eroica/route.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 12 fields
## geometry type:  LINESTRING
## dimension:      XYZ
## bbox:           xmin: 11.3464 ymin: 43.03679 xmax: 11.61873 ymax: 43.49132
## z_range:        zmin: 129.4 zmax: 620
## geographic CRS: WGS 84

r <- st_coordinates(route)[, 1:2] %>% 
  Line() %>% 
  Lines("Eroica")
r <- SpatialLines(list(r))
crs(r) <- crs(route)
plot_area(dem, route, void = TRUE)
scalebar(d = 10000, 
         type = "bar", 
         below = "kilometers",
         label = as.character(c(0, 5000, 10000) / 1000),
         div = 4)

I force the DEM on the route only:

el <- mask(dem, r) %>% 
  crop(extent(route))

The energy landscape for a cyclist of 60 kg cycling at an overall speed of 30 km/h can then be computed using to use the model from (Prampero et al. 1979) by specifying method = "cycling":

en <- enerscape(el, 60, "kcal", 4, method = "cycling", v = 30)

##  - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)

##   | Calculating slope

##   | Calculating work

##   | Calculating conductance (1 / work)

##  - Do not use slope with negative values for distance calculations

The total energy costs of cycling L’Eroica are then derived by summing the energy expenditures for the whole route:

costs <- sum(values(mask(en$rasters$Work, el)), na.rm = TRUE)
costs

## [1] 5782.049

To get the energy costs for cyclists of different weights and for other cycling speed, I make a for loop, e.g.:

library(foreach)
res <- foreach (weight = seq(50, 100, by = 10),
                .combine = "rbind") %do% {
                  res <- foreach (speed = seq(10, 50, by = 5),
                                  .combine = "rbind") %dopar% {
                                    en <- enerscape(el, 
                                                    weight, 
                                                    "kcal", 
                                                    4, 
                                                    method = "cycling", 
                                                    v = speed)
                                    cbind(speed,
                                          sum(values(mask(en$rasters$Work, el)),
                                              na.rm = TRUE))
                                  }
                  cbind(weight, res)
                } %>% 
  as_tibble()

I can then plot the energy costs for different body weight and cycling speed:

res %>% 
  mutate(weight = as.factor(weight)) %>% 
  ggplot() +
  geom_line(aes(speed, V3, col = weight)) +
  geom_violin(aes(speed, V3, group = as.factor(speed)), alpha = 0, draw_quantiles = 0.5) +
  geom_point(aes(speed, V3, col = weight), size = 1) +
  xlab("Speed (km / h)") +
  ylab("Energy consumed (kcal)") +
  scale_color_brewer("Weight (kg)", palette = "RdGy") +
  theme_classic() +
  theme(panel.grid.major.y = element_line(color = "grey80", size = 0.25))

and the difference in minimum and maximum costs at same speed, but different weight:

res %>% 
  group_by(speed) %>% 
  summarize(Min = min(V3),
            Max = max(V3)) %>% 
  mutate(Diff = Max - Min) %>% 
  ggplot() +
  aes(speed, Diff) +
  geom_smooth(alpha = 0.2, col = "tomato") +
  xlab("Speed (km / h)") +
  ylab("Maximum energy costs\ndifference (kcal)") +
  geom_point() +
  theme_classic()

## `summarise()` ungrouping output (override with `.groups` argument)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'